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Abstract 

In  turbulence  applications,  strongly  imposed  no-slip  conditions  often  lead  to  inaccurate 
mean  flow  quantities  for  coarse  boundary-layer  meshes.  To  circumvent  this  shortcoming, 
weakly  imposed  Dirichlet  boundary  conditions  for  fluid  dynamics  were  recently  introduced 
in  [4].  In  the  present  work,  we  propose  a  modification  of  the  original  weak  boundary  condi¬ 
tion  formulation  that  consistently  incorporates  the  well-known  “law  of  the  wall”.  To  com¬ 
pare  the  different  methods,  we  conduct  numerical  experiments  for  turbulent  channel  flow 
at  Reynolds  number  395  and  950.  In  the  limit  of  vanishing  mesh  size  in  the  wall-normal 
direction,  the  weak  boundary  condition  acts  like  a  strong  boundary  condition.  Accordingly, 
strong  and  weak  boundary  conditions  give  essentially  identical  results  on  meshes  that  are 
stretched  to  better  capture  boundary  layers.  However,  on  uniform  meshes  that  are  incapable 
of  resolving  boundary  layers,  weakly  imposed  boundary  conditions  deliver  significantly 
more  accurate  mean  flow  quantities  than  their  strong  counterparts.  Hence,  weakly  imposed 
boundary  conditions  present  a  robust  technique  for  flows  of  industrial  interest,  where  op¬ 
timal  mesh  design  is  usually  not  feasible  and  resolving  boundary  layers  is  prohibitively 
expensive.  Our  numerical  results  show  that  the  formulation  that  incorporates  the  law  of  the 
wall  yields  an  improvement  over  the  original  method. 
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1  Introduction 


In  computational  fluid  dynamics  formulations  that  employ  continuous  representa¬ 
tion  of  the  fields,  Diriehlet  boundary  eonditions  are  typieally  imposed  by  speeifying 
the  nodal  values  of  the  solution.  This  amounts  to  so-ealled  “strong  satisfaetion”  of 
the  boundary  eonditions.  In  flow  eomputations,  strongly  imposed  no-slip  eonditions 
often  lead  to  inaeeurate  mean  flow  quantities  for  insuffieiently  fine  boundary-layer 
meshes.  Reeently,  Bazilevs  and  Hughes  [4]  proposed  to  satisfy  Diriehlet  bound¬ 
ary  conditions  in  a  weak  sense  rather  than  strongly.  To  this  end,  the  variational 
equations  are  augmented  by  terms  that  enforce  the  Diriehlet  conditions  weakly  as 
Euler-Lagrange  eonditions.  Thus,  the  funetions  representing  the  diserete  solution 
are  not  required  to  satisfy  the  Diriehlet  conditions  explicitly.  It  was  found  that  for 
the  linear  advection-diffusion  equation  it  is  precisely  the  weak  Diriehlet  bound¬ 
ary  conditions  that  are  able  to  mitigate  or  even  entirely  eliminate  oseillations  due 
to  unresolved  boundary  layers  as  well  as  to  improve  the  aeeuraey  in  the  regions 
away  from  the  layers.  Moreover,  numerieal  results  for  low  Reynolds  number  flows 
eomputed  on  eoarse  meshes  demonstrated  that  weak  no-slip  boundary  eonditions 
provide  a  signifieant  inerease  in  aeeuraey  over  their  strong  eounterparts. 

In  the  present  work,  we  revisit  the  weak  Diriehlet  condition  formulation.  Although 
the  design  of  the  boundary  condition  is  based  on  numerical  rather  than  physieal 
considerations,  the  weak  treatment  seems  to  behave  like  a  wall  function.  To  exploit 
this  link  with  wall  modeling,  we  propose  a  modifieation  of  the  original  formula¬ 
tion  that  consistently  ineorporates  the  well-known  “law  of  the  wall”,  an  empirical 
relation  between  the  near-wall  fluid  veloeity  and  the  distanee  from  the  wall  that 
is  eommonly  assumed  to  hold  for  a  broad  range  of  Reynolds  numbers  [21].  We 
eombine  the  weakly  imposed  boundary  condition  formulation  with  residual-based 
turbulenee  modeling,  whieh  is  a  new  paradigm  for  eomputing  turbulent  flows  in- 
trodueed  in  [6,  15]  and  further  developed  in  [1].  To  eompare  the  different  Diriehlet 
boundary  eondition  formulations,  we  assess  their  performanee  on  turbulent  ehan- 
nel  flows  at  medium-to-high  Reynolds  numbers.  These  numerical  test  eases  are 
more  challenging  than  the  ones  eonsidered  previously  in  [4]  due  to  the  increased 
Reynolds  number.  In  the  limit  of  vanishing  mesh  size  in  the  wall-normal  direetion, 
the  weak  formulation  aets  like  a  strong  formulation.  Aeeordingly,  strong  and  weak 
formulations  give  essentially  identieal  results  on  stretched  meshes  that  are  designed 
to  better  resolve  the  boundary  layer.  However,  on  meshes  that  are  uniform  also  in 
the  wall-normal  direetion,  weakly  imposed  Diriehlet  boundary  conditions  deliver 
signifieantly  more  aceurate  mean  flow  quantities  than  their  strong  eounterparts. 
This  faet  makes  the  weakly  enforeed  boundary  eondition  formulations  attraetive 
for  computing  flows  of  industrial  interest,  allowing  one  to  avoid  the  eostly  resolu¬ 
tion  of  boundary  layers  without  compromising  the  aeeuraey  of  large-seale  features. 
We  also  find  that  the  weak  formulation  modified  to  ineorporate  the  law  of  the  wall 
provides  an  improvement  over  the  original  formulation.  Throughout  this  work,  the 
spatial  diseretization  makes  use  of  the  Isogeometric  Analysis  approach  [2,  3,  8, 16]. 
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The  paper  is  organized  as  follows.  In  Seetion  2,  we  deseribe  the  weak  formula¬ 
tion  of  the  eontinuous  problem  for  the  ineompressible  Navier-Stokes  equations. 
We  then  state  the  diserete,  residual-based  variational  multiseale  formulation  of  the 
problem  with  no-slip  Diriehlet  boundary  eonditions  imposed  weakly.  In  Seetion  3, 
we  deseribe  the  new  formulation  with  weakly  imposed  boundary  eonditions  that 
ineorporates  the  law  of  the  wall  by  appropriately  modifying  the  boundary  terms  of 
the  original  weak  boundary  eondition  formulation.  In  Seetion  4,  we  show  numeri- 
eal  results  for  an  equilibrium  turbulent  ehannel  flow  at  Reynolds  numbers  395  and 
950  based  on  frietion  veloeity.  In  all  oases,  we  use  meshes  for  our  eomputations 
that  are  orders  of  magnitude  ooarser  than  the  ones  employed  in  high-fidelity  Direet 
Numerioal  Simulation  (DNS);  see  [9,  20].  In  Seetion  5,  we  draw  oonolusions. 


2  Weak  Imposition  of  Diriehlet  Boundary  Conditions  for  Incompressible  Navier- 
Stokes  Equations 


2. 1  Continuous  problem 


We  begin  by  eonsidering  a  weak  formulation  of  the  Ineompressible  Navier-Stokes 
equations.  Let  V  denote  the  trial  solution  and  weighting  funetion  spaees,  whieh  are 
assumed  to  be  the  same.  We  also  assume  u  =  0  on  T  and  p{t)  dfl  =  0  for  all 
t  G  ]0,  T[.  The  variational  formulation  is  stated  as  follows:  Find  U  =  {u,p}  G  V 
sueh  that  VVE  =  {w,  q}  G  V, 


B{W,U)  =  {W,F) 


(1) 


where 


B{W,U) 


( Vm,  u  0  u)^  +  (g,  V  ■  u)n  -  ( V  ■  m,  p)^  (2) 
+  2z/V*u)^ , 


and 


{W,F)  =  {w,f)n.  (3) 

Variational  equations  (l)-(3)  imply  satisfaetion  of  the  linear  momentum  equations 
and  of  the  ineompressibility  eonstraint,  namely 

C{u,p)-f  =  0  in  fl,  (4) 

V  •  u  =  0  in  fl,  (5) 
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where 


C{u,p)  =  —  +  V  ■  (it  0  n)  +  Vp  —  V  ■  (2z/V®ii). 


(6) 


We  also  introduee  the  “advective”  form  of  the  above  operator 

Ou 

Cadv{u,p)  =  —  +  It  ■  Vn  +  Vp  -  Z/An, 


(7) 


whieh  is  obtained  from  (6)  by  using  the  ineompressibility  eonstraint  in  the  adveetive 
term  and  in  the  viseous  stress  term. 


2.2  Discrete  formulation 

Below,  we  reeall  the  discrete  variational  formulation  of  the  incompressible  Navier- 
Stokes  equations  with  weakly  imposed  Dirichlet  boundary  conditions;  see  also  [4]. 

Let  be  decomposed  into  rtez  elements,  which  induces  the  decomposition  of  V  into 
rieft  boundary  faces.  We  approximate  (l)-(3)  by  the  following  variational  problem 
over  the  finite  element  spaces:  Find  =  {ii^,p^}  G  V^,  ix^  ■  n  =  0  on  F  such 
that  \/W^  =  {lu^,  G  lu^  ■  n  =  0  on  F, 


(8) 


+  ■  Vm"  +  yq^}TM,Cadv{u\p^)  -  f)n^ 


e=l 

riel 


+  ■  i^W^T}TM,Cadv{u’",P^)  -  f)ne 


e=l 

n-el 


®rM  {Cadv{u^,P^)  -  /})Oe 


e=l 


e=l 


-  ■  n)r,nr 


b=l 

'^eb 


-  ■n,u^  -  0)r,nr 


b=l 


with  the  following  definitions 


tm  := 


(9) 
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and 


T-c  :=  {g  ■  TmO)  ^ , 

where  G  is  a  seeond-rank  metrie  tensor 

\ox  ox 


(10) 


(11) 


gi  is  a  veetor  obtained  by  summing  G  on  its  first  index  as 


S=(9).  =  E(G)ji. 


(12) 


and  ^  is  the  inverse  Jaeobian  of  the  element  mapping  between  the  parent  and  the 
physieal  domain.  In  (8)-(9),  hb  is  the  wall-normal  element  mesh-size,  and  Cj,  Ct 
and  Cl  are  positive  constants. 

Remarks 


(1)  The  above  formulation  makes  use  of  a  Residual-Based  Multiscale  Method 
(see  e.g.  [1,6,  15]),  which  is  based  on  the  Variational  Multiscale  Formulation 
(see  e.g.  [10-14,  17]).  These  residual-based  methods  possess  a  dual  nature: 
on  the  one  hand  they  are  bona-fide  LES-like  turbulence  models,  and  on  the 
other  hand  they  may  be  thought  of  as  stabilized  methods,  such  as  SUPG  [5], 
extended  to  the  nonlinear  realm. 

(2)  The  last  three  terms  of  (8)  pertain  to  the  weak  enforcement  of  the  no-slip 
condition,  as  presented  in  [4] .  We  choose  to  enforce  the  normal  component  of 
the  no-slip  boundary  condition,  that  is,  the  no-penetration  condition,  strongly 
on  the  trial  and  weighting  function  spaces. 

(3)  In  the  case  of  strongly  imposed  no-slip  conditions,  the  last  three  terms  of  (8) 
vanish. 


3  Weakly  Imposed  No-Slip  Dirichlet  Boundary  Conditions  Based  on  a  Wall 
Function  Formulation 


In  this  section,  we  revisit  weakly  imposed  Dirichlet  boundary  conditions  and  pro¬ 
pose  a  modification  of  the  original  formulation  presented  in  the  previous  section. 
This  modification  draws  on  the  knowledge  of  the  fluid  behaviour  in  the  vicinity  of 
the  wall  in  the  regime  of  fully  developed  turbulence.  In  what  follows,  we  reformu¬ 
late  the  weakly  imposed  Dirichlet  condition  in  a  way  that  is  consistent  with  the  idea 
of  wall  modeling. 
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In  engineering  praetiee  it  is  often  of  interest  to  aeeurately  resolve  large-seale  flow 
features  rather  than  fine-seale  eomponents.  It  is  typieally  not  the  detailed  features 
of  the  boundary  layer  turbulenee  that  are  relevant  for  the  applieation,  but  their  ef- 
feet  on  the  overall  flow  behavior.  This  faet  ean  be  aeeounted  for  by  wall  modeling, 
in  whieh  the  no-slip  Diriehlet  boundary  eondition  is  replaeed  by  a  traetion  Neu¬ 
mann  boundary  eondition;  see  for  example  [19,  p.  47].  In  the  direetion  tangent  to 
the  wall  a  shear  stress  is  speeified  by  adding  the  following  term  to  the  variational 
formulation 


^eb 

E 

6=1 


u 


\u 


h\  I /Tbnr  ) 


(13) 


where  1 1  ■  1 1  denotes  Euelidean  length.  The  magnitude  of  the  wall  shear  stress 
is  eonsistent  with  the  so-ealled  law  of  the  wall.  This  “law”  is  an  empirieal  relation 
between  the  mean  fluid  speed  and  the  normal  distanee  to  the  wall.  Among  the  many 
available  parameterizations  we  employ  the  one  given  by  Spalding  [21] 


2  6  ;  ’ 


(14a) 


where  and  denote  the  distanee  from  the  wall  and  the  mean  fluid  speed,  re- 
speetively,  expressed  in  non-dimensional  wall  units 


+  yu* 

y  ■=  — , 

V 


(14b) 


u  '  :  = 


(14c) 


In  Eqs.  (14b)  and  (14c),  u*  is  the  friction  velocity,  y  is  the  vertical  distance  to  the 
wall,  is  the  velocity  parallel  to  the  wall,  and  x  =  0.4  and  B  =  5.5. 


Upon  rearranging  terms  in  (13),  and  dropping  the  sum  over  the  element  boundaries 
for  brevity,  the  “penalty”  structure  of  (13)  becomes  apparent,  that  is 


acting  as  a  penalty  parameter.  Based  on  this  observation,  we  propose  to  mod¬ 
ify  the  original  weak  boundary  condition  formulation  (8)  as  follows:  Eind  = 
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{it^,  p^}  G  V^,  u^  n  =  0  on  r  such  that  yw’^  =  {lu^,  q’^}  G  V^,  w^  n  =  0  on  F, 


(17) 


»lei 


+  ■  Vw^  +  Vg'^jrM,  -  /)o, 

e=l 
riel 

+  £({'“^  ■  -  f)n, 

e=l 
riel 

-Y,{'^W\  TM{CaiM,  P”)  -  /}  0  '^M  {^adv{_'^  •)  P 

e=l 
riel 

+  Y,(^  ■  w\tcV  ■  u%, 

e=l 
^eb 

-  2vV"v!^  ■  n)r,nr 

b=i 
l^eb 

-  ^(2i/V®iu^  ■n,u^  -  0)r,nr 


-m 


b=l 


=  0. 


Variational  equation  (17)  differs  from  (8)  only  in  the  last  term  on  the  left-hand  side, 
and  it  may  be  thought  of  as  a  generalization  of  (8).  Seleeting  ?/  to  be  a  eonstant 
multiple  of  the  wall-normal  mesh  size  hb,  that  is,  y  =  hb/Cl,  and  letting  hb  go 
to  zero,  the  Spalding  equation  (14a)  reduees  to  y+  =  M+,  whieh  is  a  well-known 
parameterization  of  the  viseous  sublayer.  In  this  limit,  tb  beeomes  independent  of 
the  slip  veloeity  and  takes  on  the  expression 


V 

Tb  =  - 

y 


(18) 


Thus,  we  reeover  the  original  weak  formulation  (8).  This,  in  turn,  implies  that  the 
formulation  (17)  inherits  all  the  attributes  of  the  original  formulation  (8)  in  this 
limit.  Conversely,  when  the  mesh  size  hb  is  large,  tb  deviates  from  (18). 


Algorithm  1  outlines  a  Newton  proeedure  to  determine  tb  from  given  u^,  hb  and 
u  in  aeeordanee  with  the  law-of-the-wall  equation  (14a).  This  proeedure  is  lo- 
eal  to  eaeh  boundary-faee  integration  point.  Expression  (18)  with  y  =  hb/Cl  is 
used  to  initialize  tb-  In  ease  u^,  hb  and  u  eorrespond  to  the  viseous  sublayer,  the 
law-of-the-wall  equation  (14a)  is  satisfied  by  the  initial  values,  and  no  iteration  is 
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necessary.  Algorithm  1  makes  use  of  the  Jacobian 


dr 


hh 


dr, 


B 


2vCi 


(19) 


1  +  Xe  —  1  — 


-3/2 


xw  - 


hi  1 1/2 


where  r  ■= 


f{u^)  is  the  residual  of  the  Spalding  equation  (14a). 
Alg.  1:  Algorithm  for  computing  r^. 


1 .  Initialize  iteration  counter:  i  =  0 

2.  Initialize  TB,i  =  Cl ^ 

3-  yt  =  ut  = 

4-  A  =  yf  -  f{ut' 

5.  While(|ri|  >  TOL)  Do 


hb 

\u^ 


--1/2  1 1  hi  1 1/2 


6. 

7. 

8. 

9. 

10. 

11. 

12. 

13.  Enddo 


Build  Jacobian:  -§^\i  according  to  (19) 

Solve  for  increment:  Ars,i-i-i  =  —  (^|i) 
Update:  rs,i+i  =  Ts^i  +  ^TB,i+i 

\u^\\^C 
h||l/2 

fiut+i) 


-1 


„.+  _  hb  CC 

//*+!  vC[ 

O 


u. 


-1/2 
i+l 

A+i  =  r/i+i 

i  =  i  +  1 


4  Numerical  experiments  for  turbulent  channel  flow 

4. 1  Problem  setup 


To  investigate  the  performance  of  the  weak  boundary  condition  formulations,  we 
conduct  numerical  experiments  for  turbulent  channel  flow  at  Reynolds  numbers 
Rct  =  395  and  Rcr  =  950,  with  Re^  based  on  the  friction  velocity  and  the  channel 
half  width.  We  compare  the  results  with  the  formulation  that  imposes  the  no-slip 
condition  strongly.  To  assess  the  accuracy  of  our  methods,  we  compare  our  results 
to  the  DNS  results  of  [20]  for  Rcr  =  395  and  [9]  for  Rcr  =  950. 

The  problem  setup  is  shown  in  Figure  1 .  The  flow  is  driven  by  a  pressure  gradient  in 
the  stream- wise  direction.  At  the  computational  domain  boundary,  periodic  bound¬ 
ary  conditions  are  imposed  in  both  stream- wise  and  span- wise  directions,  whereas  a 
homogeneous  Dirichlet  boundary  condition  is  applied  in  the  wall-normal  direction. 
Stream-wise  and  span-wise  directions  are  commonly  referred  to  as  homogeneous 
directions. 
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Flow  driven  by  pressure  gradient 


Fig.  1.  Turbulent  channel  flow.  Problem  setup. 

The  spatial  discretization  is  comprised  of  quadratic  spline  functions  that  are  C^- 
continuous  at  knots.  This  type  of  discretization  is  commonly  employed  in  isoge¬ 
ometric  analysis  [16].  The  semi-discrete  equations  are  advanced  in  time  using  the 
generalized-a  method  [7,  18].  We  employ  meshes  that  are  uniform  in  all  directions. 
For  comparison  we  also  use  meshes  that  are  stretched  in  the  wall-normal  direction 
to  cluster  points  near  the  boundary  layer.  The  stretching  is  obtained  by  distributing 
the  knots  according  to  a  hyperbolic  tangent  function  such  that  the  first  knot  lies  at 
1.3,  which  is  typical  of  Large  Eddy  Simulation  (LES)  computations.  Details 
of  the  computational  setup  are  shown  in  Table  1.  Moreover,  we  set  Ct  =  4,  C/  =  36 
and  Cl  =  4. 


Table  1 

Details  of  the  computational  setup.  Lx,y,2  denotes  the  length  of  the  channel  in  the  stream- 
wise,  wall-normal  and  span-wise  direction,  Nf,i  is  the  number  of  elements  in  the  domain, 
^x,y,z  is  the  number  of  basis  functions  in  the  stream-wise,  wall-normal  and  span-wise 
direction,  fx  is  the  forcing  in  the  stream-wise  direction,  and  u  denotes  kinematic  viscosity. 


Lx 

Ly 

L, 

Nel 

Nx 

Ny 

N, 

fx 

V 

Re  =  395 

27r 

2 

323 

32 

34 

32 

3.372040  •  10-3 

1.47200  •  10-"^ 

Re  =  950 

dvr 

2 

643 

64 

66 

64 

2.630991  •  10-3 

0.53992  •  lO-'^ 

Numerical  results  for  all  cases  are  reported  in  the  form  of  statistics  of  the  mean 
velocity  and  root-mean-square  of  the  velocity  fluctuations.  Statistics  are  computed 
by  sampling  the  velocity  field  at  the  mesh  knots  and  averaging  the  solution  in  time 
as  well  as  in  the  stream-wise  and  span-wise  directions.  The  mean  velocity  is  typi¬ 
cally  referred  to  as  the  primary  statistic,  while  the  fluctuations  are  called  secondary 
statistics.  It  is  generally  acknowledged  that  accuracy  of  the  fluctuations  is  more 
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difficult  to  achieve  than  aeeuraey  of  the  mean  veloeity.  Results  are  presented  in 
non-dimensional  wall  units. 


4.2  Turbulent  channel  flow  at  Re^  =  395 


Our  eomputations  are  earried  out  on  a  mesh  of  32^  elements.  This  diseretization 
gives  32  basis  funetions  in  the  homogeneous  direetions  and  34  basis  funetions  in  the 
wall-normal  direetion  due  to  the  open  knot  veetor  eonstruetion  (see  [16]  for  details). 
In  terms  of  the  number  of  degrees-of-freedom,  this  type  of  resolution  is  typieal  of 
LES  at  Reynolds  number  395.  The  domain  size  is  27r,  2,  and  2/37r  in  the  stream- 
wise,  wall-normal,  and  span-wise  direetions,  respeetively.  The  eorresponding  DNS 
eomputation  was  earried  out  on  a  domain  of  the  same  size  and  the  diseretization 
used  256  x  193  x  192  speetral  funetions  in  the  stream-wise,  wall-normal  and  span- 
wise  direetion,  respeetively. 

Figures  2  and  3  show  statisties  of  the  eomputations  on  the  stretehed  and  uniform 
meshes,  respeetively. 

On  the  stretehed  mesh,  both  the  mean  flow  and  the  fluetuations  are  in  very  good 
agreement  with  the  DNS  (see  Figure  2).  In  faet,  the  aeeuraey  of  the  results  is  vir¬ 
tually  that  of  a  speetral  eomputation,  although  simple  quadratie  spline  funetions 
with  loeal  support  are  used  instead  of  speetral  basis  funetions.  Results  obtained 
with  strongly  and  weakly  imposed  no-slip  eonditions  praotieally  eoineide,  whieh  is 
eonsistent  with  the  faet  that  the  weak  boundary  eondition  formulation  reduees  to 
the  strong  one  in  the  limit  of  vanishing  mesh  size.  The  newly  proposed  formulation 
that  ineorporates  wall  modeling  gives  slightly  more  aeeurate  stream-wise  veloeity 
fluetuations  than  the  other  formulations. 

In  eontrast,  on  the  uniform  mesh,  the  methods  perform  differently  (see  Figure 
3).  Plaeing  the  first  knot  at  23,  we  intentionally  saerifice  the  resolution  of 

the  boundary  layer.  This  leads  to  a  gross  overestimation  of  the  mean  flow  for  the 
strongly  enforeed  Dirichlet  boundary  eondition  formulation.  Note  that,  on  the  other 
hand,  the  mean  veloeity  for  both  weak  formulations  agrees  very  well  with  the  DNS 
result.  This  shows  that  weak  boundary  eonditions  are  eapable  of  alleviating  the 
gross  inaeeuraey  indueed  by  insuffieient  near- wall  resolution.  This  superior  robust¬ 
ness  despite  “poor”  mesh  design  makes  the  new  method  attraetive  for  industrial 
applieations.  We  also  note  that  the  wall  funetion  formulation  is  slightly  more  aeeu¬ 
rate  for  mean  flow  veloeity  than  the  original  weak  boundary  eondition  formulation. 
Despite  the  large  differenee  in  the  mean  flow,  the  seeondary  statisties  in  the  eore  of 
the  ehannel  for  the  uniform  mesh  eases  are  very  similar  for  all  formulations  eonsid- 
ered.  In  the  near  wall  region,  the  differenees  in  the  fluetuations  obtained  with  the 
various  methods  are  more  pronounced. 

Comparing  the  results  obtained  on  stretched  and  on  uniform  meshes,  we  observe 
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Fig.  2.  Turbulent  channel  flow  at  Rcr  =  395  computed  on  a  stretched  mesh.  Top:  Mean 
stream-wise  velocity  plotted  versus  wall  distance  in  wall  units.  Bottom:  Stream-wise, 
wall-normal  and  span-wise  velocity  fluctuations  plotted  versus  wall  distance  in  wall  units. 
Formulation  with  no-slip  boundary  conditions  enforced  strongly  (□),  weakly  according  to 
original  methodology  (8)  (A),  and  weakly  based  on  the  wall  function  (17)  (-i-). 

that  the  secondary  statistics  for  the  uniform  mesh  simulations  are  not  quite  as  ac¬ 
curate  as  those  for  the  stretched  grid  case,  although  the  quality  of  the  results  is 
still  good.  One  may  thus  conclude  that  in  the  core  of  the  channel  the  effect  of  the 
mesh  design  on  the  fluctuations  is  more  pronounced  than  the  effect  of  the  boundary 
conditions. 
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Fig.  3.  Turbulent  channel  flow  at  Rcr  =  395  computed  on  a  uniform  mesh.  Top:  Mean 
stream-wise  velocity  plotted  versus  wall  distance  in  wall  units.  Bottom:  Stream-wise, 
wall-normal  and  span-wise  velocity  fluctuations  plotted  versus  wall  distance  in  wall  units. 
Formulation  with  no-slip  boundary  conditions  enforced  strongly  (□),  weakly  according  to 
original  methodology  (8)  (A),  and  weakly  based  on  the  wall  function  (17)  (-i-). 


Figure  4  shows  the  stream- wise  veloeity  eontours  at  an  instant  in  time,  computed 
on  a  uniform  mesh  with  weak  boundary  conditions  employing  the  wall  function 
formulation.  Note  the  presence  of  velocity  fluctuations  of  considerable  magnitude 
at  the  “no-slip”  wall  (top  surface  of  the  box). 
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Fig.  4.  Turbulent  channel  flow  at  Rer  =  395.  Snapshot  of  stream-wise  velocity  contours. 
4.3  Turbulent  channel  flow  at  Rcr  =  950 


For  the  computations  at  Rct  =  950,  a  mesh  of  64^  elements  is  used  with  64 
basis  functions  in  the  homogeneous  directions  and  66  basis  function  in  the  wall- 
normal  direction  due  to  an  open  knot  vector  construction.  The  domain  size  is  Ttt,  2, 
and  4/37r  in  the  stream- wise,  wall-normal  and  span- wise  directions,  respectively. 
The  corresponding  DNS  used  a  domain  size  of  Stt  x  2  x  Svr  with  a  resolution  of 
3072  X  385  x  2304  spectral  functions  in  the  stream-wise,  wall-normal  and  span- 
wise  directions.  Note  that  our  resolution  per  unit  domain  length  is  a  factor  of  about 
25  coarser  in  the  stream-wise  direction,  a  factor  of  6  coarser  in  the  wall-normal 
direction,  and  a  factor  of  24  coarser  in  the  span-wise  direction.  Hence,  the  adopted 
discretization  is  significantly  coarser  than  what  is  typically  used  for  an  LES-type 
computation. 

Figures  5  and  6  show  statistics  of  the  computations  on  stretched  and  uniform  meshes, 
respectively. 

On  a  stretched  mesh,  the  differences  between  the  weak  and  the  strong  boundary- 
condition  formulations  are  negligible  due  to  the  small  near- wall  mesh  size  in  the 
wall-normal  direction  (see  Figure  5).  All  methods  fail  to  accurately  represent  the 
mean  flow  velocity.  Moreover,  the  stream- wise  velocity  fluctuations  are  inaccurate 
in  the  near- wall  region  but  are  quite  accurate  in  the  core  of  the  channel.  The  ve¬ 
locity  fluctuations  in  the  remaining  directions  are  in  very  good  agreement  with  the 
DNS.  The  good  agreement  of  the  velocity  fluctuations  with  the  DNS  despite  the 
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Fig.  5.  Turbulent  channel  flow  at  Rer  =  950  computed  on  a  stretched  mesh.  Top:  Mean 
stream-wise  velocity  plotted  versus  wall  distance  in  wall  units.  Bottom:  Stream-wise, 
wall-normal  and  span-wise  velocity  fluctuations  plotted  versus  wall  distance  in  wall  units. 
Formulation  with  no-slip  boundary  conditions  enforced  strongly  (□),  weakly  according  to 
original  methodology  (8)  (A),  and  weakly  based  on  the  wall  function  (17)  (-i-). 

discrepancy  in  the  mean  flow  velocity  is  somewhat  surprising.  The  inability  to  ac¬ 
curately  capture  the  mean  velocity  illustrates  the  limitations  of  the  strong  boundary- 
condition  method  for  high  Reynolds  number  wall-bounded  flows. 

On  the  uniform  mesh,  with  the  first  knot  at  t/+  30,  the  strong  boundary  condition 

formulation  gives  an  even  greater  over-prediction  of  the  mean  velocity  than  for 
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Fig.  6.  Turbulent  channel  flow  at  Rcr  =  950  computed  on  a  uniform  mesh.  Top:  Mean 
stream-wise  velocity  plotted  versus  wall  distance  in  wall  units.  Bottom:  Stream-wise, 
wall-normal  and  span-wise  velocity  fluctuations  plotted  versus  wall  distance  in  wall  units. 
Formulation  with  no-slip  boundary  conditions  enforced  strongly  (□),  weakly  according  to 
original  methodology  (8)  (A),  and  weakly  based  on  the  wall  function  (17)  (-i-). 


the  stretched  mesh  computations;  compare  Figures  5  and  6.  In  contrast,  both  weak 
boundary-condition  formulations  deliver  a  result  of  reasonable  accuracy  for  such 
a  coarse  discretization,  similar  to  the  case  of  Rcr  =  395.  The  mean  flow  is  only 
slightly  over-predicted  when  compared  to  the  DNS.  Also  note  that  the  uniform 
discretization  does  not  exhibit  as  severe  an  overshoot  in  the  stream- wise  velocity 
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fluctuations  near  the  wall  as  for  the  stretched  mesh.  However,  away  from  the  wall, 
the  fluetuations  are  slightly  less  accurate  on  the  uniform  mesh  than  on  the  stretehed 
mesh.  This  is  eonsistent  with  the  results  for  the  Re^  =  395  case. 

Figure  7  shows  the  stream-wise  veloeity  eontours  at  an  instant  in  time,  eomputed 
on  a  uniform  mesh  with  weak  boundary  eonditions  employing  the  wall  funetion 
formulation.  Note  the  presenee  of  veloeity  fluctuations  of  considerable  magnitude 
at  the  “no-slip”  wall  (top  surfaee  of  the  box).  Also  note  that  the  turbulent  struetures 
for  the  Rer  =  950  channel  are  more  fine-grained  than  the  ones  for  the  Rcr  =  395 
ehannel  due  to  the  increased  Reynolds  number  (eompare  Figures  4  and  7). 


5  Conclusions 


In  this  work,  we  proposed  a  new  variational  formulation  of  the  incompressible 
Navier-Stokes  equations  that  enforees  Diriehlet  no-slip  boundary  eonditions  weakly. 
Motivated  by  the  observation  that  weak  imposition  of  Diriehlet  boundary  condi¬ 
tions  generally  seems  to  behave  like  a  wall  funetion,  the  proposed  formulation  is 
based  on  the  so-ealled  law  of  the  wall.  We  eombined  the  weak  boundary  eondi- 
tion  formulation  with  residual-based  turbulence  modeling.  We  compared  the  per¬ 
formance  of  the  different  boundary  condition  formulations  based  on  numerical  re¬ 
sults  for  turbulent  ehannel  flow  at  medium-to-high  Reynolds  number.  We  found 
that  the  weakly  imposed  boundary  condition  formulation  that  incorporates  the  law 
of  the  wall  provides  an  improvement  over  the  original  weak  boundary  eondition 
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formulation.  In  the  limit  of  vanishing  mesh  size  in  the  wall-normal  direetion,  both 
weak  boundary  eondition  formulations  aet  like  a  strong  formulation.  Aeeordingly, 
weak  and  strong  boundary  eondition  formulations  give  essentially  identieal  results 
on  stretehed  meshes  that  are  designed  to  better  eapture  the  boundary  layer.  How¬ 
ever,  on  eoarse,  uniform  meshes,  weakly  imposed  boundary  eonditions  deliver  a 
signifieantly  more  aeeurate  mean  flow  veloeity  than  their  strong  eounterpart.  In 
this  respeet,  the  eombination  of  residual-based  turbulenee  modeling  and  weak  im¬ 
position  of  the  no-slip  eondition  aets  like  a  RANS-type  model  in  the  sense  that 
it  produees  aeeurate  mean  flow  quantities  on  meshes  that  are  too  eoarse  for  eon- 
ventional  LES  simulations.  This  result  makes  weakly  imposed  Diriehlet  boundary 
eondition  formulations  attraetive  for  eomputing  flows  of  industrial  interest,  avoid¬ 
ing  the  eostly  resolution  of  boundary  layers  without  eompromising  the  aeeuraey  of 
large-seale  flow  features.  Given  that  the  weak  boundary  eondition  formulation  be¬ 
haves  like  its  strong  eounterpart  on  fine  meshes  and  delivers  superior  aeeuraey  on 
eoarse  meshes,  this  suggests  the  use  of  this  method  as  a  general  strategy  for  enfore- 
ing  wall  boundary  eonditions  in  finite-element  flow  eomputations.  The  additional 
eost  due  to  weak  enforeement  of  the  boundary  eonditions  is  negligible  beeause  the 
eorresponding  integrals  are  evaluated  only  over  the  Diriehlet  portion  of  the  domain 
boundary. 

Regarding  the  role  of  weak  versus  strong  boundary  eondition  formulation  in  the 
eontext  of  residual-based  turbulenee  modeling,  our  results  for  the  Rcr  =  395  ehan- 
nel  flow  showed  that  the  residual  based  formulation  with  standard  strongly  im¬ 
posed  boundary  eonditions  gives  remarkably  aeeurate  results  for  a  well  designed, 
stretehed  mesh  with  a  resolution  that  is  typieal  of  LES.  However,  for  flows  at  higher 
Reynolds  number,  sueh  as  Rcr  =  950  eomputed  on  a  stretehed  mesh  with  a  resolu¬ 
tion  eorresponding  to  eoarse  LES  /  fine  RANS,  the  aeeuraey  of  the  residual-based 
approaeh  deteriorates.  The  laek  of  aeeuraey  in  predieting  the  mean  flow  veloeity 
derives  from  the  inability  to  resolve  the  turbulent-boundary-layer  flow  struetures 
on  a  mesh  that  is  insuffieiently  fine.  Our  results  demonstrated  that  this  diffieulty 
ean  be  elegantly  eireumvented  by  eombining  the  residual-based  formulation  with  a 
weakly  enforeed  no-slip  boundary  eondition. 
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